{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Combustion Model and Ignition Delay Demo\n",
    "\n",
    "Written by A. Mark Payne and Alon Grinberg Dana for presentation at ICCK 2017"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## User input:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "fuel = 'OC'\n",
    "equivalence_ratio = 1.0\n",
    "\n",
    "temperature = 1500.0  # (K)\n",
    "pressure = 1.0        # (atm)\n",
    "sim_time = 2          # (ms)\n",
    "top_sens = 10         # number of top sensitive reactions and thermo to display\n",
    "\n",
    "rmgpy_path = '../rmg.py' # Change to your rmg.py path\n",
    "\n",
    "from IPython.display import display, Image\n",
    "from rmgpy.molecule import Molecule\n",
    "\n",
    "fuel_molecule = Molecule(smiles=fuel)\n",
    "print(\"The fuel molecule is:\")\n",
    "display(fuel_molecule)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## RMG's input file:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "fuel_molecule = Molecule(smiles=fuel)\n",
    "nC = int(fuel_molecule.get_num_atoms('C'))\n",
    "nH = int(fuel_molecule.get_num_atoms('H'))\n",
    "nO = int(fuel_molecule.get_num_atoms('O'))\n",
    "\n",
    "fuel_stoich = equivalence_ratio/(nC+(nH/4.0)-(nO/2.0))\n",
    "\n",
    "input_file = f'''\n",
    "# Data sources\n",
    "database(\n",
    "    thermoLibraries = ['BurkeH2O2','primaryThermoLibrary','thermo_DFT_CCSDTF12_BAC','DFT_QCI_thermo','FFCM1(-)'],\n",
    "    reactionLibraries = ['BurkeH2O2inN2','FFCM1(-)'],\n",
    "    seedMechanisms = [],\n",
    "    kineticsDepositories = ['training'],\n",
    "    kineticsFamilies = 'default',\n",
    "    kineticsEstimator = 'rate rules',\n",
    ")\n",
    "\n",
    "# List of species\n",
    "species(\n",
    "    label='fuel',\n",
    "    reactive=True,\n",
    "    structure=SMILES('{fuel}'),\n",
    ")\n",
    "\n",
    "species(\n",
    "    label='O2',\n",
    "    reactive=True,\n",
    "    structure=SMILES('[O][O]'),\n",
    ")\n",
    "\n",
    "species(\n",
    "    label='N2',\n",
    "    reactive=False,\n",
    "    structure=SMILES('N#N'),\n",
    ")\n",
    "\n",
    "species(\n",
    "    label='OH',\n",
    "    reactive=True,\n",
    "    structure=SMILES('[OH]'),\n",
    ")\n",
    "\n",
    "# Reaction system\n",
    "simpleReactor(\n",
    "    temperature=({temperature!s},'K'),\n",
    "    pressure=({pressure!s},'atm'),\n",
    "    initialMoleFractions={{\n",
    "        'fuel': {fuel_stoich!s},\n",
    "        'O2': 1,\n",
    "        'N2': 3.76,\n",
    "    }},\n",
    "    terminationTime=({sim_time/1000.0},'s'),\n",
    "    sensitivity=['OH'],\n",
    "    sensitivityThreshold=0.01,\n",
    ")\n",
    "\n",
    "simulator(\n",
    "    atol=1e-16,\n",
    "    rtol=1e-8,\n",
    "    sens_atol=1e-6,\n",
    "    sens_rtol=1e-4,\n",
    ")\n",
    "\n",
    "model(\n",
    "    toleranceKeepInEdge=0,\n",
    "    toleranceMoveToCore=0.1,\n",
    "    toleranceInterruptSimulation=0.1,\n",
    "    maximumEdgeSpecies=100000,\n",
    "    filterReactions=True,\n",
    "    maxNumObjsPerIter=2,\n",
    "    terminateAtMaxObjects=True,\n",
    "    maxNumSpecies=50,\n",
    ")\n",
    "\n",
    "#pressureDependence(\n",
    "#        method='modified strong collision',\n",
    "#        maximumGrainSize=(0.5,'kcal/mol'),\n",
    "#        minimumNumberOfGrains=250,\n",
    "#        temperatures=(298,2500,'K',10),\n",
    "#        pressures=(0.5,3,'bar',5),\n",
    "#        interpolation=('Chebyshev', 6, 4),\n",
    "#        maximumAtoms=16,\n",
    "#)\n",
    "\n",
    "options(\n",
    "    units='si',\n",
    "    generateOutputHTML=True,\n",
    "    generatePlots=False,\n",
    "    saveEdgeSpecies=False,\n",
    "    saveSimulationProfiles=True,\n",
    ")\n",
    "\n",
    "generatedSpeciesConstraints(\n",
    "    allowed=['input species','seed mechanisms','reaction libraries'],\n",
    "    maximumCarbonAtoms=5,\n",
    "    maximumOxygenAtoms=2,\n",
    "    maximumNitrogenAtoms=0,\n",
    "    maximumSiliconAtoms=0,\n",
    "    maximumSulfurAtoms=0,\n",
    "    maximumHeavyAtoms=6,\n",
    "    maximumRadicalElectrons=2,\n",
    "    allowSingletO2=False,\n",
    ")\n",
    "\n",
    "\n",
    "'''\n",
    "\n",
    "import os\n",
    "import shutil\n",
    "directory = './rmg_demo'\n",
    "if os.path.exists(directory):\n",
    "    shutil.rmtree(directory)\n",
    "os.mkdir(directory)\n",
    "input_path = os.path.join(directory, 'input.py')\n",
    "with open(input_path,'w') as f:\n",
    "    f.write(input_file)\n",
    "print('Created RMG input file at ' + os.path.abspath(input_path))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run RMG:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "import datetime\n",
    "import subprocess\n",
    "start = time.time()\n",
    "\n",
    "# Execute RMG job\n",
    "subprocess.check_call(['python', rmgpy_path, input_path])\n",
    "\n",
    "end = time.time()\n",
    "print('Total simulation time: ' + str(datetime.timedelta(seconds=round(end-start))))\n",
    "\n",
    "with open(os.path.join(directory, 'RMG.log'),'r') as f:\n",
    "    begin = False\n",
    "    for line in f:\n",
    "        if 'MODEL GENERATION COMPLETED' in line:\n",
    "            begin = True\n",
    "        if begin:\n",
    "            print(line.strip())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Run the generated model (using RMG's Cantera functions):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "from rmgpy.chemkin import load_chemkin_file\n",
    "from rmgpy.tools.canteraModel import Cantera, get_rmg_species_from_user_species\n",
    "from rmgpy.species import Species\n",
    "import time\n",
    "\n",
    "chem_path = os.path.join(directory, 'chemkin')\n",
    "species_list, reaction_list = load_chemkin_file(os.path.join(chem_path, 'chem_annotated.inp'),\n",
    "                                              os.path.join(chem_path, 'species_dictionary.txt'),\n",
    "                                              os.path.join(chem_path, 'tran.dat'))\n",
    "\n",
    "fuel_species = Species(smiles=fuel)\n",
    "O2_species = Species(smiles='[O][O]')\n",
    "N2_species = Species(smiles='N#N')\n",
    "OH_species = Species(smiles='[OH]')\n",
    "species_dict = get_rmg_species_from_user_species([fuel_species, O2_species, N2_species, OH_species], species_list)\n",
    "\n",
    "reactor_type_list = ['IdealGasReactor']\n",
    "reaction_time_list = ([sim_time], 'ms')\n",
    "\n",
    "mol_frac_list=[{species_dict[fuel_species]: fuel_stoich,\n",
    "                species_dict[O2_species]: 1,\n",
    "                species_dict[N2_species]: 3.76}]\n",
    "T_list = ([temperature],'K')\n",
    "P_list = ([pressure],'atm')\n",
    "\n",
    "job = Cantera(species_list=species_list, reaction_list=reaction_list, output_directory=directory)\n",
    "job.load_chemkin_model(os.path.join(chem_path, 'chem_annotated.inp'), os.path.join(chem_path, 'tran.dat'))\n",
    "job.generate_conditions(reactor_type_list, reaction_time_list, mol_frac_list, T_list, P_list)\n",
    "\n",
    "alldata = job.simulate()\n",
    "print(\"Simulation Completed\")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Plot:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "############### Settings ###############\n",
    "fsize = (8,4) # Change to make the figure fit on your screen\n",
    "########################################\n",
    "\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "from rmgpy.tools import plot as rmg_plot\n",
    "from operator import itemgetter\n",
    "%matplotlib notebook\n",
    "\n",
    "times = alldata[0][0].data\n",
    "temperatures = alldata[0][1][0].data\n",
    "pressures = alldata[0][1][1].data\n",
    "\n",
    "dpdt = (pressures[1:] - pressures[:-1]) / (times[1:] - times[:-1])\n",
    "idi = next(i for i,d in enumerate(dpdt) if d==max(dpdt))\n",
    "ign_delay_time = times[idi]\n",
    "\n",
    "for spc in range(len(alldata[0][1][:])):\n",
    "    if alldata[0][1][spc].label == str(species_dict[fuel_species]):\n",
    "        Fuel_idx = spc\n",
    "    if alldata[0][1][spc].label == str(species_dict[OH_species]):\n",
    "        OH_idx = spc\n",
    "\n",
    "for i in range(len(alldata[0][1][Fuel_idx].data)):\n",
    "    if alldata[0][1][Fuel_idx].data[i]<0.001:\n",
    "        Fuel_Depletion_Time = times[i]\n",
    "        break\n",
    "\n",
    "files = os.listdir(os.path.join(directory, 'solver'))\n",
    "\n",
    "sensitivity_file = [f for f in files if ('sensitivity' in f) and ('.csv' in f)][0]\n",
    "SA_time, SA_data = rmg_plot.parse_csv_data(os.path.join(directory, 'solver', sensitivity_file))\n",
    "\n",
    "time_error = 1\n",
    "\n",
    "for i in range(len(SA_time.data)):\n",
    "    if abs(SA_time.data[i]-ign_delay_time)<time_error:\n",
    "        ign_delay_idx = i\n",
    "        time_error = abs(SA_time.data[i]-ign_delay_time)\n",
    "\n",
    "Gidx = 0\n",
    "for i in range(len(SA_data[:])):\n",
    "    if \"G\" in SA_data[i].label:\n",
    "        if not Gidx:\n",
    "            Gidx = i\n",
    "        SA_data[i].label = SA_data[i].label.split('G')[1][1:-1]\n",
    "    else:\n",
    "        SA_data[i].label = SA_data[i].label.split(' ')[1]\n",
    "\n",
    "rank1 = []\n",
    "for n in range(Gidx):\n",
    "    rank1.append(abs(SA_data[n].data[ign_delay_idx]))               # list of max SA range for each rxn\n",
    "num1 = np.linspace(0,len(rank1)-1,len(rank1))        # will be used to get the order of reactions by rank1\n",
    "num1 = zip(rank1,num1)\n",
    "num1 = sorted(num1, key=itemgetter(0),reverse=True)\n",
    "SA_k_data = []\n",
    "SA_k_label = []\n",
    "for i in range(min(top_sens, Gidx)):\n",
    "    SA_k_data.append(SA_data[int(num1[i][1])].data[ign_delay_idx])     # make sorted lists size topSA of SA values and rxns labels\n",
    "    SA_k_label.append(SA_data[int(num1[i][1])].label)\n",
    "        \n",
    "rank2 = []\n",
    "for n in range(len(SA_data)-Gidx):\n",
    "    rank2.append(abs(SA_data[n+Gidx].data[ign_delay_idx]))     # list of max SA range for each rxn\n",
    "num2 = np.linspace(0,len(rank2)-1,len(rank2))        # will be used to get the order of reactions by rank1\n",
    "num2 = zip(rank2,num2)\n",
    "num2 = sorted(num2, key=itemgetter(0),reverse=True)\n",
    "SA_G_data = []\n",
    "SA_G_label = []\n",
    "for i in range(min(top_sens, len(SA_data)-Gidx)):\n",
    "    SA_G_data.append(SA_data[int(num2[i][1])+Gidx].data[ign_delay_idx])     # make sorted lists size topSA of SA values and rxns labels\n",
    "    SA_G_label.append(SA_data[int(num2[i][1])+Gidx].label)\n",
    "    \n",
    "print(\"Ignition delay time is {0:.4f} ms\".format(ign_delay_time * 1000))\n",
    "        \n",
    "plt.rcParams['axes.labelsize'] = 18\n",
    "plt.rcParams['xtick.labelsize'] = 12\n",
    "plt.rcParams['ytick.labelsize'] = 12\n",
    "plt.rcParams['figure.autolayout'] = True\n",
    "\n",
    "plt.style.use('ggplot')\n",
    "plt.style.use('seaborn-pastel')\n",
    "\n",
    "fig = plt.figure(figsize=fsize)\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(alldata[0][0].data*1000, alldata[0][1][Fuel_idx].data,'-o')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('$Y_{Fuel}$')\n",
    "plt.title('Fuel profile')\n",
    "plt.xlim([0,2000*ign_delay_time])\n",
    "\n",
    "max_oh = max(alldata[0][1][OH_idx].data)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(alldata[0][0].data*1000, alldata[0][1][OH_idx].data,'-o')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('$Y_{OH}$')\n",
    "plt.title('OH profile')\n",
    "plt.xlim([0,2000*ign_delay_time])\n",
    "plt.arrow(0, alldata[0][1][OH_idx].data[idi], ign_delay_time*1000, 0, width=max_oh*0.01, head_width=max_oh*0.05, head_length=ign_delay_time*120, length_includes_head=True, color='r', shape='full')\n",
    "plt.annotate(r'$Ignition Delay: \\tau_{ign}$', xy=(0,0), xytext=(0, alldata[0][1][OH_idx].data[idi]+0.0005), fontsize=10);\n",
    "\n",
    "fig = plt.figure(figsize=fsize)\n",
    "\n",
    "plt.subplot(1,2,1)\n",
    "plt.plot(alldata[0][0].data*1000, temperatures,'-o')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Temperature (K)')\n",
    "plt.title('Temperature')\n",
    "plt.xlim([0,2000*ign_delay_time])\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.plot(alldata[0][0].data*1000, pressures,'-o')\n",
    "plt.xlabel('Time (ms)')\n",
    "plt.ylabel('Pressure (Pa)')\n",
    "plt.title('Pressure')\n",
    "plt.xlim([0,2000*ign_delay_time])\n",
    "\n",
    "fig = plt.figure(figsize=fsize)\n",
    "plt.barh(np.arange(min(top_sens, Gidx)), SA_k_data, 1/1.5, color=\"blue\")\n",
    "plt.gca().invert_yaxis()\n",
    "plt.xlabel(r'Sensitivity: $\\frac{\\partial\\:\\ln{[OH]}}{\\partial\\:\\ln{k}}$');\n",
    "plt.rcParams.update({'axes.labelsize': 20})\n",
    "plt.yticks(np.arange(min(top_sens, Gidx)),SA_k_label)\n",
    "plt.title(\"[OH] sensitivity to kinetics\")\n",
    "\n",
    "fig = plt.figure(figsize=fsize)\n",
    "plt.barh(np.arange(min(top_sens, len(SA_data)-Gidx)), SA_G_data, 1/1.5, color=\"blue\")\n",
    "plt.gca().invert_yaxis()\n",
    "plt.xlabel(r'Sensitivity: $\\frac{\\partial\\:\\ln{[OH]}}{\\partial\\:G_i}$ $[mol/kcal]$');\n",
    "plt.rcParams.update({'axes.labelsize': 20})\n",
    "plt.yticks(np.arange(min(top_sens, len(SA_data)-Gidx)),SA_G_label)\n",
    "plt.title(\"[OH] sensitivity to thermo\")\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cantera as ct\n",
    "gas = ct.Solution(os.path.join(directory, 'cantera', 'chem.cti'))\n",
    "comp = str(species_dict[fuel_species])+\":\"+str(fuel_stoich)+\",\"+str(species_dict[O2_species])+\":1,\"+str(species_dict[N2_species])+\":3.76\"\n",
    "gas.TPX = temperature, pressure, comp\n",
    "reactor = ct.IdealGasConstPressureReactor(gas)\n",
    "network = ct.ReactorNet([reactor])\n",
    "network.advance(ign_delay_time)\n",
    "ROP_C = ct.ReactionPathDiagram(gas, 'C')\n",
    "\n",
    "from PIL import Image as PILimg\n",
    "ROP1 = plt.subplot(1,1,1)\n",
    "dot_file = os.path.join(directory, 'cantera', 'rxnpathC.dot')\n",
    "img_file = os.path.join(directory, 'cantera', 'rxnpathC.png')\n",
    "ROP_C.title = 'Reaction path diagram following C'\n",
    "ROP_C.threshold = 0.01\n",
    "ROP_C.label_threshold = 0.01\n",
    "ROP_C.show_details = True\n",
    "ROP_C.write_dot(dot_file)                                              # write dot file\n",
    "os.system('dot {0} -Tpng -o{1} -Gdpi=300'.format(dot_file, img_file))  # write png file\n",
    "fullpath = os.getcwd() + '/' + img_file\n",
    "display(Image(fullpath))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:rmg_env]",
   "language": "python",
   "name": "conda-env-rmg_env-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.4"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
